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Abstract 

To  determine  the  potential  impacts  of  fuel  cells  on  future  distribution  system,  dynamic  models  of  fuel  cells  should  be  created,  reduced 
in  order,  and  scattered  throughout  test  feeders. 

This  paper  presents  the  implementation  of  an  efficient  method  for  computing  low-order  linear  system  models  of  solid  oxide  fuel  cells 
(SOFCs)  from  time  domain  simulations.  The  method  is  the  Box-Jenkins  algorithm  for  calculating  the  transfer  function  of  a  linear  system 
from  samples  of  its  input  and  output. 

©  2003  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Exhibiting  the  dynamic  influences  of  solid  oxide  fuel  cell 
(SOFC)  on  the  distribution  grid  requires  the  use  of  a  large 
dynamic  model  [1].  Since  SOFCs  will  be  proliferated,  it  is 
necessary  to  reduce  the  model  order  of  each  SOFC  system 
to  enable  computational  analysis. 

The  computation  of  linear  system  models  of  power  sys¬ 
tems  from  time  domain  simulations  is  a  topic  of  considerable 
practical  interest.  This  interest  is  motivated  by  the  insight 
into  the  dynamic  interactions  among  power  system  com¬ 
ponents  that  can  be  obtained  from  a  linear  representation. 
Linear  models  allow  for  the  application  of  linear  analysis 
techniques  to  complement  the  information  obtained  from 
nonlinear  time  domain  simulations  and  often  allow  for  a 
better  understanding  of  the  system  dynamic  characteristics 
than  that  obtained  from  the  inspection  of  time  simulations 
alone.  Although  the  nonlinear  nature  of  a  SOFC  must  be 
recognized,  in  many  cases  a  linearized  system  representa¬ 
tion  allows  for  a  more  efficient  means  of  analysis. 

Several  techniques  for  computing  state  space  matrices 
and  transfer  function  realizations  of  power  systems  from 
time  domain  data  have  been  proposed  in  recent  years. 
These  techniques  include  the  Prony  method  which  is  based 
on  fitting  a  weighted  sum  of  exponential  terms  to  a  given 
signal  [2],  methods  based  on  FFT  analyses  [3],  and  the 
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eigensy stem  realization  algorithm  [4] .  In  addition,  methods 
designed  to  estimate  the  frequency  response  of  a  system 
have  also  been  proposed  [5].  The  Prony  method  is  perhaps 
the  method  for  which  more  extensive  results  and  applica¬ 
tions  have  been  documented.  Its  application  to  the  analysis 
and  control  of  electromechanical  oscillations  has  shown  the 
value  of  deriving  linear  models  from  time  domain  simula¬ 
tions  and  measured  data  [6]. 

This  paper  presents  the  application  of  the  Autoregres¬ 
sion  with  exogenous  signal  (ARX)  identification  algorithm 
to  compute  low-order  system  models,  suitable  for  analysis 
and  control  design  [7-9].  This  algorithm  consists  of  a  sim¬ 
ple  procedure  for  calculating  the  transfer  function  of  a  linear 
system  from  samples  of  its  input  and  output. 

Using  MATLAB/Simulink  [10],  a  dynamic  model  of  a 
SOFC-penetrated  distribution  system  is  created. 

This  paper  is  structured  as  follows.  Section  2  presents  a  re¬ 
view  of  the  SOFC.  Section  3  introduces  the  utility-connected 
inverter  control.  Some  basic  concepts  of  ARX  models  are 
described  in  Section  4.  Section  5  compares  the  response  of 
identified  system  versus  the  response  of  the  actual  system. 
Section  6  depicts  some  simulation  results.  Finally,  conclu¬ 
sions  are  presented  in  Section  7. 


2.  Solid  oxide  fuel  cell 

Most  likely,  fuel  cell  will  be  a  dominant  distributed  en¬ 
ergy  resource.  The  SOFCs  are  dynamic  devices  and  when 
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Nomenclature 

Fuel  cell 

Eo 

ideal  standard  potential 

F 

Faraday’s  constant 

he 

fuel  cell  current 

anode  valve  constant 

Kr2 

valve  molar  constant  for  hydrogen 

Kr20 

valve  molar  constant  for  water 

Kq2 

valve  molar  constant  for  oxygen 

Kr 

constant  (=Ar0/4JF) 

^h2o 

molecular  mass  of  hydrogen 

«h2 

number  of  hydrogen  moles  in  the 

No 

anode  channel 

number  of  cells  in  series  in  the  stack 

Pi 

partial  pressure 

P 

real  power 

p* 

set  point  for  the  real  power 

4 

input  fuel  flow 

4 

output  fuel  flow 

<?h2 

fuel  flow  that  reacts 

r 

ohmic  loss 

m-o 

ratio  of  hydrogen  to  oxygen 

R 

universal  gas  constant 

T 

absolute  temperature 

Te 

electrical  response  time 

Tf 

fuel  processor  response  time 

u 

fuel  utilization 

Van 

volume  of  the  anode 

Vfc 

fuel  cell  voltage 

th2 

response  time  for  hydrogen  flow 

^h2o 

response  time  for  water  flow 

^o2 

response  time  for  oxygen  flow 

Inverter 

E 

load  bus  voltage 

E* 

set  point  for  the  load  bus  voltage 

Lp 

inductance 

Q 

reactive  power 

Q* 

set  point  for  the  reactive  power 

V 

inverter  output  voltage  space  vector 

XT 

reactance  ( =Lpco ) 

8p 

angle  between  i/r v  and  \/fe 

s; 

angle  reference 

te 

flux  vector  associated  with  E 

tv 

flux  vector  associated  with  V 

t*v 

flux  vector  reference 

connected  to  the  distribution  system  they  will  affect  its  dy¬ 
namic  behavior.  Hence,  researchers  have  developed  dynamic 
models  for  these  components  [11-15]. 

The  chemical  response  in  the  fuel  processor  is  usually 
slow.  It  is  associated  with  the  time  to  change  the  chemical 
reaction  parameters  after  a  change  in  the  flow  of  reactants. 


This  dynamic  response  function  is  modeled  as  a  first-order 
transfer  function  with  a  5  s  time  constant. 

The  electrical  response  time  in  the  fuel  cells  is  gener¬ 
ally  fast  and  mainly  associated  with  the  speed  at  which  the 
chemical  reaction  is  capable  of  restoring  the  charge  that  has 
been  drained  by  the  load.  This  dynamic  response  function 
is  also  modeled  as  a  first-order  transfer  function  but  with  a 
0.8  s  time  constant. 

With  aid  of  the  inverter,  the  fuel  cell  system  can  supply 
not  only  real  power  but  also  reactive  power.  Usually,  power 
factor  can  be  in  the  range  of  0. 8-1.0.  The  SOFC  system 
dynamic  model  is  given  in  Fig.  1. 

The  fuel  utilization  (U)  is  the  ratio  between  the  fuel  flow 
that  reacts  and  the  input  fuel  flow  as  follows: 


C  -  ?h5 


Typically,  an  80-90%  fuel  utilization  is  used. 

Every  individual  gas  will  be  considered  separately,  and 
the  perfect  gas  equation  will  be  applied  to  it.  Hydrogen  will 
be  considered  as  an  example 


PU2  Van  =  nu2RT 


It  is  possible  to  isolate  the  pressure  and  to  take  the  time 
derivative  of  the  previous  expression,  obtaining 


dpn2  RT 
d t  ~  VanqH2 


There  are  three  relevant  contributions  to  the  hydrogen  molar 
flow:  the  input  flow,  the  flow  that  takes  part  in  the  reaction 
and  the  output  flow,  thus 


According  to  the  basic  electrochemical  relationships,  the 
molar  flow  of  hydrogen  that  reacts  can  be  calculated  as 


Np  I 
IF 


=  2AVfc 


Returning  to  the  calculation  of  the  hydrogen  partial  pressure, 
it  is  possible  to  write 


If  it  could  be  considered  that  the  molar  flow  of  any  gas 
through  the  valve  is  proportional  to  its  partial  pressure  inside 
the  channel,  according  to  the  expressions: 


<?h2 

Pn2 


Replacing  the  output  flow  by  (7),  taking  the  Laplace  trans¬ 
form  of  both  sides  and  isolating  the  hydrogen  partial  pres¬ 
sure,  yields  the  following  expression: 


P  h2  = 


V*H2  ,  in 

TTw  ^ 


where  th2  =  Vm/ Ku2RT. 
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Fig.  1.  SOFC  system  dynamic  model. 


A  similar  operation  can  be  made  for  all  the  reactants  and 
products.  Applying  Nernst’s  equation  and  Ohm’s  law  (to 
consider  ohmic  losses),  the  stack  output  voltage  is  repre¬ 
sented  by  the  following  expression: 


V  =  Nq 


In 


n  ,-,0.5 
PKlPOi 

Pu2o 


3.  Utility-connected  inverter  control 

Control  of  the  flux  vector  has  been  shown  to  have  good 
dynamic  and  steady-state  performance  [16].  It  also  pro¬ 
vides  a  convenient  means  to  define  the  power  angle  since 
the  inverter  voltage  vector  switches  position  in  the  d-q 
plane,  whereas  there  is  no  discontinuity  in  the  inverter  flux 
vector. 

For  a  six-pulse  voltage  source  inverter  (VSI),  the  inverter 
output  voltage  space  vector  can  take  any  of  the  seven  po¬ 
sitions  in  the  plane  specified  by  the  d-q  coordinates.  The 
time  integral  of  the  inverter  output  voltage  space  vector  is 
called  the  “inverter  flux  vector”  for  short.  The  d-  and  q- axes 
components  of  the  inverter  flux  vector  xj/v  are  defined  as 

Fjdr,  \// dq=  f  Vq  dr  (10) 

00  J — 00 


The  magnitude  and  the  angle  of  xlsv  with  respect  to  the  g-axis 
are  determined  as 

IW=V^i,  <5„  =  tan”1  (11) 

The  d-  and  q- axes  components  of  the  ac  system  voltage  flux 
vector  xjfe ,  its  magnitude,  and  angle  are  defined  in  a  similar 
manner.  The  angle  between  x//v  and  x//e  is  defined  as 

8p  =  8V  —  8e  (12) 

It  is  useful  to  develop  the  power  transfer  relationships  in 
terms  of  the  flux  vectors.  The  basic  real  power  transfer  rela¬ 
tionship  for  the  control  system  of  Fig.  2  in  the  d-q  reference 
frame  is 

P  =  j(eqiq  +  edid)  (13) 

In  (13),  eq  and  ed  are  the  q-  and  d- axes  components,  respec¬ 
tively,  of  the  ac  system  voltage  vector  E.  In  addition,  iq  and 
id  are  the  components  of  the  current  vector  I.  When  iq  and 
id  are  expressed  in  terms  of  the  fluxes,  taking  into  account 
the  spatial  relationships  between  the  two  flux  vectors  and 
assuming  the  ac  system  voltage  to  be  sinusoidal,  (13)  can 
be  expressed  as 

3 

P  =  ——(o\lfe\lrvsmSp 

ZLj1 


(14) 
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Fig.  2.  Control  system  for  the  inverter. 


In  this  expression  \jfe  and  \/rv  are  the  magnitudes  of  the  ac 
system  and  the  inverter  flux  vectors,  respectively,  and  8p  is 
the  spatial  angle  between  the  two  flux  vectors,  i  is  the  fre¬ 
quency  of  rotation  of  the  two  flux  vectors.  The  expression 
for  reactive  power  transfer  can  be  derived  in  a  similar  man¬ 
ner.  This  is 

Q  =  cos  8p  -  ft]  (15) 

Eqs.  (14)  and  (15)  indicate  that  P  can  be  controlled  by  con¬ 
trolling  8P,  which  can  be  defined  as  the  power  angle,  and  Q 
can  be  controlled  by  controlling  ^rv. 

The  two  variables  that  are  controlled  directly  by  the  in¬ 
verter  are  v  and  8p.  The  vector  xj/y  is  controlled  to  have  a 
specified  magnitude  and  a  specified  position  relative  to  the 
ac  system  flux  vector  \fre. 

The  errors  between  actual  and  desired  amounts  activate 
the  remainder  of  the  firing  scheme  only  if  they  exceed  a 
threshold  value.  If  the  error  is  larger  than  the  hysteresis  band 
(whose  widths  are  A 8p  and  A\/rv)  then  a  decision  towards  a 
new  switching  sequence  is  made.  If  the  errors  are  within  their 
hysteresis  band,  the  switches  will  hold  their  current  status. 
Therefore,  the  SOFC  plant  has  two  major  control  loops: 

1.  Power  control :  done  by  adjusting  the  set  point  P*  of  the 
inverter  for  fast  transient  variations  and  fuel  flow  input 
control  for  slow  variations. 


2.  Voltage  control :  done  by  adjusting  the  set  point  E *  of 
the  inverter,  which  effects  the  magnitude  of  the  converter 
output  voltage. 

4.  Identification  algorithms 

4.1.  ARX  models 

The  most  used  model  structure  is  the  simple  linear  differ¬ 
ence  equation 

y(t)  +  a\y(t  -  1)  4 - h  anay(t  -  no) 

=  b\u(t  —  nk)  +  •  •  •  +  bnt>u(t  —  nk  —  nb  +  1)  +  eft) 

(16) 

which  relates  the  current  output  y(t)  to  a  finite  number  of 
past  outputs  y(t  —  k )  and  inputs  u(t  —  k). 

The  structure  is  thus  entirely  defined  by  the  three  integers 
na,  nb,  and  nk.  na  is  equal  to  the  number  of  poles  and  nb  —  1 
is  the  number  of  zeros,  while  nk  is  the  pure  time-delay  (the 
dead-time)  in  the  system.  For  a  system  under  sampled-data 
control,  typically  nk  is  equal  to  1  if  there  is  no  dead¬ 
time. 

For  multi-input  systems  nb  and  nk  are  row  vectors,  where 
the  ith  element  gives  the  order/delay  associated  with  the  ith 
input. 
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Dynamic  response  of  SOFC  plant 


There  are  two  methods  to  estimate  the  coefficients  a  and 
in  the  ARX  model  structure: 

Least  squares'.  Minimizes  the  sum  of  squares  of  the 
right-hand  side  minus  the  left-hand  side  of  the  expression 
above,  with  respect  to  a  and  b. 

Instrumental  variables’.  Determines  a  and  b  so  that  the 
error  between  the  right-  and  left-hand  sides  becomes 


uncorrelated  with  certain  linear  combinations  of  the 
inputs. 

4.2.  ARMAX,  output- error  and  Box-Jenkins  models 

There  are  several  elaborations  of  the  basic  ARX  model, 
where  different  disturbance  models  are  introduced.  These 


Bode  diagram 


Fig.  4.  Transfer  function  magnitude  and  phase  comparison.  Actual  system  and  identified  system. 
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include  well  known  model  types,  such  as  ARM  AX, 
output-error,  and  Box-Jenkins  [17-19]. 

A  general  input-output  linear  model  for  a  single-output 
system  with  input  u  and  output  y  can  be  written  as 


nu 


Mq)y(t )  =  y  fBi(q)Fj(q)]ui(t  -  nk, )  + 


i=  1 


C(q) 

L  £>(<?)  J 


e(t) 


(17) 


Here  wz-  denotes  input  #/,  and  A,  Bi,  C,  D ,  and  iq  are  poly¬ 
nomials  in  the  shift  operator  (z  or  q).  It  is  just  a  compact 
way  of  writing  difference  equations. 

The  general  structure  is  defined  by  giving  the  time-delays 
nk  and  the  orders  of  these  polynomials  (i.e.,  the  number  of 
poles  and  zeros  of  the  dynamic  model  from  u  to  y,  as  well 
as  of  the  disturbance  model  from  e  to  y). 

Most  often  the  choices  are  confined  to  one  of  the  following 
special  cases: 


O  DG2 


Fig.  5.  One  line  diagram  of  IEEE  13  node  feeder  with  fuel  cells. 


ARX  :  A(q)y(t)  =  B(q)u(t  —  nk)  +  e(t) 

ARM  AX  :  A(q)y(t)  =  B(q)u(t  —  nk)  +  C(q)e(t) 

u(t  —  nk)  +  eft) 


Output-error  :  y(t)  = 


B(q) 

L  m  I 


(18) 

(19) 

(20) 


Box-Jenkins  : 


~B(qy 

_  m  _ 


u{t  —  nk)  + 


(21) 


Note  that  A(q)  corresponds  to  poles  that  are  common 
between  the  dynamic  model  and  the  disturbance  model. 
Likewise  Ffq)  determines  the  poles  that  are  unique  for  the 
dynamics  from  input  #/,  and  D(q)  the  poles  that  are  unique 
for  the  disturbances. 

Although  each  method  has  a  somewhat  different  set  of 
parameters  that  a  system  analyst  can  adjust,  one  requirement 
an  identified  system  must  meet  is  that  its  response  to  a  given 
input  should  match  the  response  of  the  actual  system.  This 
practical  criterion  was  used  as  a  guideline  to  adjust  the  order 
of  the  identified  systems. 


Frequency  step  transient 


Fig.  6.  Fuel  cell  response  to  a  frequency  step  transient  at  node  634  for  the  IEEE  13  node  feeder.  Real  power. 
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Frequency  step  transient 


5.  Performance 

SOFC  models  used  in  the  distribution  system  analysis 
were  constructed  as  shown  in  the  following: 

•  There  is  one  4.16kV/480  V  transformer. 

•  All  SOFCs  were  connected  at  the  end  of  their  respective 
feeders  at  480  V. 


In  this  paper  SOFC  modeled  is  denoted  as  “the  actual 
system”.  Once  an  identified  system  was  obtained,  its  time  do¬ 
main  response  and  transfer  function  were  compared  against 
the  corresponding  quantities  of  the  actual  system.  To  ac¬ 
complish  this,  the  actual  system  was  linearized  around  an 
operating  point. 

The  results  presented  here  correspond  to  a  0.02  p.u.  by 
0.1  s  probing  pulse  into  the  real  power  block  of  SOFC  in 
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Fig.  8.  One  line  diagram  of  IEEE  34  node  feeder  with  fuel  cells. 
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Fig.  1.  The  sampling  time  was  0.01  s  and  600  points  were 
used  to  perform  the  system  identification. 

Assume  a  SOFC  is  operating  with  constant  rated  voltage 
and  power  demand  0.6  p.u.  There  is  0.3  p.u.  of  step  increase 
in  the  total  load  at  t  =  10  s. 

Fig.  3  compares  the  time  response  of  identified  system  ver¬ 
sus  the  response  of  the  actual  system.  The  identified  system 
was  obtained  using  Box-Jenkins  algorithm,  and  is  of  fourth 
order.  This  method  estimates  parameters  of  the  Box-Jenkins 
model  structure  using  a  prediction  error  method.  The  order 
of  the  identified  system  is  the  minimum  order  required  to 
obtain  a  good  time  domain  match. 

Fig.  4  compares  the  magnitude  and  phase  of  the  trans¬ 
fer  function  (s)/P(s)  of  the  identified  and  the  linearized 
actual  systems.  These  plots  show  a  very  good  match  in  the 
frequency  range. 

6.  Simulation  results 

The  IEEE  test  feeders  [20]  are  used  as  the  test  system 
to  investigate  the  dynamic  characteristics  of  the  distribution 
system  with  fuel  cells.  Figs.  5,  8  and  10  show  the  test  systems 
with  the  fuel  cells. 

In  this  paper,  all  loads  are  balanced,  and  characterized  by 
constant  power. 

The  majority  of  data  for  the  fuel  cell  model  has  been  ex¬ 
tracted  from  Kuipers  and  Singhal  [21,22],  and  a  commercial 
leaflet  describing  a  SOFC  100  kW  plant. 


Here,  all  fuel  cells  in  the  test  feeders  have  the  same  dy¬ 
namic  response  and  share  the  generation  equally.  The  IEEE 
13  node  test  feeder  is  very  small,  short  and  relatively  highly 
loaded  for  a  4.16kV  feeder.  Penetration  means  the  pro¬ 
portion  of  the  distribution  feeder  load  being  supplied  by 
SOFCs  associated  with  the  distribution  feeder  [23].  In  this 
model,  an  initial  load  of  Px  is  assumed  and  the  penetration  is 
thus, 

P 

Penetration  =  -  (22) 

P  +  Pi 

in  this  paper,  the  penetration  level  of  the  IEEE  test  feeders 
is  set  at  10%. 

The  first  controlled  transient  was  a  0.1  p.u.  step  in  fre¬ 
quency  at  the  point  of  connection  of  the  distribution  system 
at  0.5  s,  while  voltage  was  held  constant  (Figs.  6  and  7).  The 
second  controlled  transient  was  a  0.1  p.u.  step  in  voltage  at 
the  point  of  connection  of  the  distribution  system  at  0.5  s, 
while  frequency  was  relatively  stable  during  each  transient 
(Figs.  9,  11  and  12). 

The  IEEE  34  node  test  feeder  is  an  actual  feeder.  The 
feeder’s  nominal  voltage  is  24.9 kV.  It  is  characterized  by 
very  long  and  lightly  loaded. 

The  IEEE  123  node  test  feeder  operates  at  a  nominal 
voltage  of  4.16kV.  It  does  provide  voltage  drop  problems 
that  must  be  solved  with  the  application  of  voltage  regulators 
and  shunt  capacitors. 

Fig.  12  shows  the  response  to  a  10%  step  of  system 
voltage. 


Voltage  step  transient 


Fig.  9.  Fuel  cell  response  to  a  voltage  step  transient  at  node  848  for  the  IEEE  34  node  feeder. 
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Fig.  10.  One  line  diagram  of  IEEE  123  node  feeder  with  fuel  cells. 


Voltage  step  transient 


Fig.  11.  Fuel  cell  response  to  a  response  for  a  voltage  step  transient  at  node  46  for  the  IEEE  123  node  feeder. 
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Voltage  step  transient 


7.  Conclusions 

The  capability  to  calculate  low-order  equivalent  linear 
systems  from  time  domain  simulations  of  SOFC  models 
using  the  ARX  algorithm  has  been  established.  After  the 
SOFC  model  was  created,  it  was  reduced  to  transfer  func¬ 
tions  using  the  ARX  algorithm;  thus,  the  transfer  function 
(reduced- order  model)  exhibited  the  same  dynamic  response 
as  the  original  SOFC  model. 

A  significant  reduction  in  the  model  order  was  achieved. 
The  time  domain  response  of  the  identified  system  matches 
the  response  of  the  actual  system. 

Therefore,  each  SOFC  reduced-order  model  influences 
the  grid  in  the  same  manner  as  SOFC  and  loads  would, 
modulating  real  and  reactive  power  in  response  to  voltage 
and  frequency  changes  on  the  grid. 
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